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Abstract 

I Many mathematical models involve input parameters, which are not precisely known. Global 

' — ^ ■ sensitivity analysis aims to identify the parameters whose uncertainty has the largest impact on 

' the variability of a quantity of interest (output of the model). One of the statistical tools used 

to quantify the influence of each input variable on the output is the Sobol sensitivity index. We 
^ ' consider the statistical estimation of this index from a finite sample of model outputs. We study 

I asymptotic and non-asymptotic properties of two estimators of Sobol indices. These properties 

, are applied to significance tests and estimation by confidence intervals. 
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1 Introduction 



Many mathematical models encountered in applied sciences involve a large number of poorly-known 
parameters as inputs. It is important for the practitioner to assess the impact of this uncertainty 
on the model output. An aspect of this assessment is sensitivity analysis, which aims to identify 
the most sensitive parameters, that is, parameters having the largest influence on the output. In 
global stochastic sensitivity analysis (see for example [14] and references therein) the input variables 
are assumed to be independent random variables. Their probability distributions account for the 
practitioner's belief about the input uncertainty. This turns the model output into a random variable, 
whose total variance can be split down into different partial variances (this is the so-called Hoeffding 
decomposition, see [17]). Each of these partial variances measures the uncertainty on the output 
induced by each input variable uncertainty. By considering the ratio of each partial variance to the 
total variance, we obtain a measure of importance for each input variable that is called the Sobol 
index or sensitivity index of the variable [15]; the most sensitive parameters can then be identified 
and ranked as the parameters with the largest Sobol indices. 

Once the Sobol indices have been defined, the question of their effective computation or estimation 
remains open. In practice, one has to estimate (in a statistical sense) those indices using a finite 
sample (of size typically in the order of hundreds of thousands) of evaluations of model outputs [3] . 
Indeed, many Monte Carlo or quasi Monte Carlo approaches have been developed by the experimental 
sciences and engineering communities. This includes the Sobol pick-freeze (SPF) scheme (see [15, 16]). 
In SPF a Sobol index is viewed as the regression coefficient between the output of the model and its 
pick-freezed replication. This replication is obtained by holding the value of the variable of interest 
(frozen variable) and by sampling the other variables (picked variables). The sampled replications 
are then combined to produce an estimator of the Sobol index. In this paper we study very deeply 
this Monte Carlo method in the general framework where one or more variables can be frozen. This 
allows to define sensitivity indices with respect to a general random input living in a probability space 
(groups of variables, random vectors, random processes...). 

In [7] , the authors have studied the asymptotic behavior of two pick-freeze estimators of a single Sobol 
index. The results in this paper can be continued in two directions. The first direction is motivated 
by the fact that in general, so as to rank input variables according to their importance, the pratic- 
tioners jointly estimate the collection of all the first-order as well as the total Sobol indices. As these 
different estimators are dependent, the asymptotic marginal distributions are not fully informative, 
and one has to characterize the joint law of the estimators. This joint law allows, for example, to 
perform significance tests and comparisons between different indices, so as to rigorously rank the 
input variables, taking into account indices estimation errors. The second direction is motivated by 
the fact that asymptotic distributions are unattainable in practice, hence, non-asymptotic tools (such 
as concentration inequalities, and Berry-Esseen-like theorems) about the distribution of the Sobol 
indices estimators should be investigated. Such results will allow conservative certification for the 
index estimates. 

This paper is organized as follows: in Section 2, we review the Sobol pick-freeze method and give 
the estimators that are studied in the paper. In Section 3, we prove a central limit theorem which 
gives the joint asymptotic distribution of any closed Sobol index [14], which in particular can be used 
to explicit the asymptotic distribution of all first-order and total index estimators. We then apply 
this central limit theorem to significance and comparison tests on Sobol indices. Sections 4 and 5 are 
dedicated to non-asymptotic studies of the distribution of a single Sobol index estimator. These two 
sections, respectively, give concentration inequalities and Berry-Esseen bounds. All our theoretical 
results are numerically illustrated on model examples. 



2 Sobol pick freeze Monte Carlo method 
2.1 Black box model and Sobol indices 

In the whole paper, we consider a non necessarily linear regression model connecting an output y G M 
to independent random input vectors Xi, . . . Xp with for i ~ 1, . . .p^ Xi belongs to some probability 
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space Xi. We denote 

y = :=/(Xi,...,Xp) (1) 

where / is a deterministic real valued measurable function defined on X = Xi x . . . Xp. We assume 
that Y is square integrable and non deterministic (VarY ^ 0). 

Let u := {ui, . . . , u^) be k subsets of Ip := {1, . . . The vector of closed Sobol indices (see [14]) is 
then 

/Var(E(y|X„i e Ml)) Var(E(r|X,,i e Uk)y 



ou 

'-'ci 



\ Var(r) Var(r) 



As pointed out and discussed in the Introduction, Sobol indices are useful quantities widely used 
in engineering and applied sciences in the context of prioritisation of influent input variables of a 
complicated computer simulation code (see for example [14], [2[) and our paper gives a rigourous 
statistical analysis of these quantities. Notice that considering the whole vector Sq^ allows estimation 
of asymptotic confidence regions and tests for joint significance (see Section 3) . 

2.2 Monte Carlo estimation of 5": Sobol pick freeze method 

For X and for any subset v of Ip we define by the vector such that X"^ ^ Xi \i i & v and X" = X[ 
\i i ^ V where X[ is an independent copy of Xi. We then set 

The next lemma [7, Lemma 1.2[ shows how to express Sq-^ in terms of covariances. This will lead to 
a natural estimator: 

Lemma 2.1. For any u d Ip, one has 

Var(E(y|X,i e m)) = Cov(y,y"). (2) 

An estimator with a close expression has been considered in [5[. 
Notation 

From now on, we will denote Var(F) by V, Cov(F, F") by C„ and Zjv the empirical mean of any 
iV-sample (Zi, . . . , Zm) of Z. 

A first estimation for S^j. In view of Lemma 2.1, we are now able to define a first natural 
estimator of Sq^ (all sums are taken for i from 1 to N) : 

These estimators have been considered in [5], where it has been showed to be practically efficient 
estimators. 

A second estimation for 5*"^. Since the observations consist in {YijY^^'- , . . . , yi"'°)(i<i<Af), a more 
precise estimation of the first and second moments can be done and we are able to define a second 
estimator of 5'"[ taking into account all the available information. Define 



k 



The second estimator is then defined as 
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This estimator (in the fc = 1 case) was first introduced by Monod in [9] and Janon et aL studied its 
asymptotic properties (CLT, efiiciency) in [7]. In [11, 10] Owen introduces new estimators for Sobol 
indices and compares numerically their performances. The delta method can also be used on these 
pick-freeze estimators to derive their asymptotic properties. 

Remark 2.2. One could use all the information available in the sample by defining the following 
estimator: 

However, our empirical studies show that this estimator has a larger variance than 

3 Joint CLT for Sobol index estimates with applications to sig- 
nificance tests 

3.1 Main results 

Theorem 3.1. Assume i/ia< E(F'') < oo. Then: 
1. 

,ci - ^Si) „^ AA, (0, r„,5) (5) 

A/— >oo 

where Tu^s = ii^u,s)i,j)i<ij<k 

_ Cov(yy"' , yy"^ ) ~ 5'^'iCov(yy"^ ,y^)- s^l Cov(yy"' ,y^) + s^l s^\Yai{Y'^) 

[i-u,S)l,j — — , ,,2 

(Var(r)) 

2. 

^ [Tn.cx - ^a) „^ A/l- (0, Tu.t) (6) 

A' — >-CX3 

where T^,t = ((ru,T)i j)i<, j<fc with 

_ Cov(yy"' , yy ) - s'^^Coy{yv^ , m") - s^^ Cov(yy"' , a/") + s^i Var(Af ") 

(ru^rjij 2 ■ 

(Var(y)) 

3.2 Some particular cases 

1. Assume k = p, u = ({1}, . . . , {p}) and E{Y'^) < oo. We denote y/-'^ by Y^^ . Here 

^ / Var(E(y|Xi)) Var(E(y|Xp)) 
V Var(y) Var(y) 



and 



Tu 
AT ni — 



The CLT becomes 



where T^^t = i(^u,T)i.j)i<ij<k '^ith 

(Var(y))^ (ru,T)ij = Cov(yy',yy^)-S'^iCov(yy^Af")-S'^iCov(yy', Af")+s'^i5^jVar(Af^ 
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2. We can obviously have a CLT for any index of order 2. Indeed if we take k ~ 1 and E 
{1, ...,p}'^ with i ^ j and u = {ij}. We get = i (y + y") and A/" = + (^")^); 
thus 

^„ ^ Var(E(y|X.,X,-)) ^^^^ ^„ ^ iEy^y^''-{^EiY^ + Yny 

Var(y) ^-^^ _^^(y2 + (yu)2)_(_i_^(y^ + yu))2- 

The CLT becomes 
with 

(Var(y))' (Fu^t) = Var(yy") - 25SiCov(yy", y^) + (Var(y2) + cov(y2, (y")')) . 

3. One can also straightforwardly deduce the joint distribution of the vector of all indices of order 
2. For example, if p = 3 take k ~ 3 and u = ({1, 2}, {1, 3}, {2, 3}) and apply Theorem 3.1. 



3.3 Proof of Theorem 3.1 

Since qj and are invariant by any centering (translation) of the li's and Y^' 's for j = 1, . . . , fc, 

we can simplify the next calculations translating by E{Y). For the sake of simplicity, Yi and y^"^ now 
denote the centered random variables. 

Proof of (5) : 

Recall that 

Let Wi = (yiy"' , j = 1, . . . , fc, y„ y"^ j = 1 . . . , fc, y/)* (i = l, . . .) and g the mapping from R^fc+s to 
M'= defined by 



g{xi,...,Xk,y,yi,-.-,yk,z) 



xi - yyi xk - yyk 
z — y' z — y' 



Let E denote the covariance matrix of Wt . The vectorial central limit theorem implies that 

We then apply the so-called Delta method [17] to W and g so that 

^{giWN)-9{nW))) 4 AA(0, Jg(E(VF))I]Jg(E(VF))*) 

with Jg(E(H/)) the Jacobian of g at point E{W). 
Define (gi, . . . , gfc) ;= For i = 1, . . . , fc, j = 1, . . . , k, 



dy 
09 j 

with 6^^^ = 1 and = if i 7^ j. Thus Fu,s = Jg(E(W^))SJg(E(M^))* is as stated in Theorem 3.1. 



Proof of (6) : 

The proof is similar to the one of (5). We now define Wi = {YiY^\j = 1, . . . , k,Yi,Y^\ j 
1 . . . , k, (y")2)*. We apply the delta method to g from M^fc+a jj^^Q jgfe defined by 

g{xi,...,xk,y,yi,.-.,yk,z) 



„ _ i y+yi+---+yk \ ■y ^ ( v+vi+---+Vk \ 

V fe+i ) ^ \ k+i ) 
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For i = j = l,...,fc, 



'-§t{nw)) = - 



^5 



V 



3.4 Significance tests 

In order to simplify the notation we will write the vectors as column vectors. In this section, we 
give a general procedure to build significance tests of level a and then illustrate this procedure on two 
examples. 

Let u := . . . , Uk) so that for any i = 1, . . . , fc, is a subset of Ip :— {1, . . . Similarly, let 
V := (ui, . . . , vi) and w := (wi, . . . , wi) be I be so that for any i = 1, . . . , C and Wi C Ip. 

Consider the following general testing problem 

Ho : S*") = and 5^1 = Sqi against Hi : Hq is not true. 
Remark 3.2. Note that one can also test 

Hq : Sqi < s against Hi : 5"; > s, 

or 

Ho '■ Sqi < Sqi against Hi : Sqi > 
Appling Theorem 3.1 we have 

G^:=Vn(( ^^'C' ) - f c^v ^"'.w ) ) ^k+i (0, r) . (7) 

Since we have an explicit expression of F we may build an estimator Fjv of F thanks to empirical 
means. Note that (F^v)^ converges a.s. to F. Define 

CU 



Gn VN 



QV CW 



Then: 



Gn — Gn — f ''^ 



QV CW 



Corollary 3.3. Under Hq, Gn 4 A/'fe+;(0,F). 
[/nder i/i, |Gjv(l)| + |GAr(2)| oo. 

This corollary allows us to construct several tests. It is a well-known fact that in the case of a vectorial 
null hypothesis "there exists no uniformly most powerful test, not even among the unbiased tests" (see 
Chapter 15 in [17]). In practice, we return to the dimension 1 introducing a function F : M'^+' — > R and 
testing HoiF) : F{h) = (respectively Hi{F) : F{h) ^ 0) instead Ho : h = Q (resp. Hi : h ^ 0). 
The choice of a reasonable test "depends on the alternatives at which we wish a high power". 

Remark 3.4. // we take as test statistic Tn = AGn where A is a linear form defined on ]R'+'^, under 
Ho,Tn Af (0, ATA'). Replacing T by Fiv CLnd using Slutsky^s lemma we get 

N—>QO 

{ATn^)-^''^Tn 4 AA(0,1). 

A^— >-oo 

Thus we reject Ho if {AT nA')^^/'^Tn > where is the 1 — a quantile of a standard Gaussian 
random variable. 

One can have a similar result when A is not anymore linear but only by applying the so-called 
Delta method. 
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3.4.1 Numerical applications: toy examples 

Example 1 In this first toy example, we compare 5 different test statistics through their power 
function. Let X = {Xi,X2) -7V(0,/2), and 

Y = f{X) = XiXi + X1X2 + X2X1X2, 

with 2X1 + X2 — I- We consider here the following testing problem 

Ha : 5*01 ^ S^^ = Xj = against Hi: Xi ^ 0. 

Then, computations lead to 

r(l, 1) = r(2, 2) = 3 - 2A? - llA^ + 24A? - 24A^ 

r(2, 1) = r(l, 2) = -7Xj + 24A? - 24Af . 

The Gaussian limit in Theorem 3.1 is A/'2(0, 3/^2) under Hq while it is asymptotically distributed as 
7V2(0,r) under Hi. 

Test 1: we take as test statistic T^ i — GAr(l) + Gn{2). 

Under Hq, Tns ^ A/'(0,6) so we reject H^ if Tat,! > Za where Zq,/-\/6 is the (1 — a) quantile of 

N^oo 

a standard Gaussian random variable. While under Hi, following the procedure of Remark 3.4 with 
A = {1 1). 

(Tjv,i-2ViVA?)/(2[r(l,l) + r(l,2)])i/2 4 AA(0,1). 

It is then easy to compute the theoretical power function. In Figure 1 we plot this theoretical function 
called true power fct tl and the empirical power function called estimated power fct tl. To 
compute the empirical power function we didn't assume the knowledge of the matrix T nor the one 
of the function /. 

Test 2: since the Sobol indices are non negative, the testing problem is naturally unilateral. However 
in view of more general contexts we introduce the test statistic Tn^2 = |G'Ar(l)| + |Gjv(2)|. We reject 
Hf) if Tn,2 > Za where Za/VS is the (1 — a) quantile of the random variable having 

^e-"'/4cI>(V^/2)lRJw) 

as density ($ being the distribution function of a standard Gaussian random variable). Under Hi, the 
power function of T/v 2 and the limit variance are estimated using Monte Carlo technics. In Figure 1 
we plot this empirical power function called estimated power fct t2. 

Test 3: in the same spirit, we introduce the test statistic Tat. 3 = |GAr(l) + GAr(2)|. We reject Hq 
if Tn^3 > Za where Za/\/6 is the (1 — a/2) quantile of a standard Gaussian random variable. Under 
Hi, the power function of T/v.s and the limit variance are estimated using Monte Carlo technics. In 
Figure 1 we plot this empirical power function called estimated power fct t3. 

Test 4: we use the L norm and consider Tna = (GAr(l))2 + (GAr(2))2. Under Hq, Tna/5 4 X2(2) 

N—yoo 

SO we reject Hq if T/v, 4 > Za where Za/S is the (1 — a) quantile of a X2 random variablewith 2 degrees 
of freedom. Under Hi, the power function of ^i^d the limit variance are estimated using Monte 
Carlo technics. We plot this empirical power function in Figure 1 called estimated power fct t4. 

Test 5: we use the infinity norm and consider Tjv_5 = TOaa;(|GAr(l)|; |GAr(2)|). We reject Hq if 
Tn.5 > Za where Za/VS is the [1 + y/l — a]/2 quantile of a standard Gaussian random variable. Under 
Hi, the power function of r/v,5 and the limit variance are estimated using Monte Carlo technics. In 
Figure 1 we plot this theoretical function called true power fct t5 and the empirical power function 
called estimated power fct t5. 

In Figure 1 we thus present the plot of the different power functions for = 100, 500 and 1000. 
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True power fct t1 








Estimated power fct t1 








1 1 1 Estimated power fct t2 








• • • Estimated power fct t3 








^1*^ ^1^ * Estimated power fct t4 








C O Estimated power fct t5 

















0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 

lambdaj, N=100 




0.0 ^ ■ ' ' ' ' ' ' ' ■ ' ' ' ' 1 

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 

lambdaj , N=500 




0.0 ^ • ' • ' • ' • ' ■ ' • ' • 1 

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 

lambdaj ,N= 1000 



Figure 1: Power functions 



Example 2 Let X = (Xi, X2, X3) - 7V(0, h), 2Af + = 1 and 

Y = f{X) = Xi{X2 + X3) + X2X1X2 




Let us test if Xi has any influence ie Hq 
Theorem 3.1 we easily get 



cfi} _ 



c{l,2} _ 



o{2} 



and S'qj'^^ — 



_ c{3} 

— ^ci 



Applying 



G 



N 



ql;2 _ C2 



q1,3 



c3 




N- 



AA3(o,r). 



Here under Hq the covariance limit F in Theorem 3.1 is the identity matrix. Under Hi we use its 
explicit expression given in Theorem 3.1 to compute an empirical estimator Fjv. We compare Test 
1, Test 3, Test 4 and Test 5 defined in the previous example. We present in Figure 2 the plot of the 
different estimated power functions for N = 100, 500 and 1000. 

Figures 1 and 2 show, as expected, that increasing leads to a steeper power function (hence, a 
better discrimination between the hypothesis), and that the estimated power function gets closer to 



8 




0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 

lambdaj, N=100 




0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 

lambdaj, N=1000 



Figure 2: Estimated power functions for different values of N. 

the true one. We also see that no test is the most powerful, uniformly in Ai, in accordance with the 
theory quoted above. 

Ishigami function The Ishigami model [6] is given by: 

Y = f{Xi,X2,X3) = sinXi + 7sin2X2 + 0.1X|sinXi (8) 

for {Xj)j^i^2,3 are i.i.d. uniform random variables in [— 7r;7r]. Exact values of these indices are 
analytically known: 

= 0.3139, 4f = 0.4424, S^f = 0. 
We perform simulations in order to show that our test procedure allows us to recover the fact that 
= 0, even for relatively small values of A^. In Table 1, we present the simulated confidence levels 
obtained for N G {10, 50, 100, 500, 1000} by the following procedure. For each value of N, we use a 
1000 sample to estimate the confidence level and we repeat this scheme 20 times. We give in Table 
1 the minimum, the mean and the maximum of these 20 distinct simulated values of the confidence 
levels. 
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N 


Mill 


Mean 


Max 


10 


0.041 


0.0463 


0.048 


50 


0.042 


0.0482 


0.050 


100 


0.044 


0.0489 


0.051 


500 


0.047 


0.0510 


0.053 


1000 


0.049 


0.0510 


0.055 



Table 1: Results for the Ishigami function 



3.4.2 Numerical applications: a real test case 

It is customary in aeronautics to model the fuel mass needed to link two fixed countries with a 
commercial aircraft by the Breguet formula: 

Mfuel = (Mempty + Mpioad) (e^^v.^"" _ 1^ . (9) 

See [13] for the description of the model with more details. 
The fixed variables are 

• Mernpty ■ Empty Weight = basic weight of the aircraft (excluding fuel and passengers) 

• Mpioad ■ Payload = maximal carrying capacity of the aircraft 

• g : Gravitational constant 

• Ra : Range = distance traveled by the aircraft 
The uncertain variables are 

• V : Cruise speed = aircraft speed between ascent and descent phase 

• F : Lift-to-drag ratio = aerodynamic coefficient 

• SFC : Specific Fuel Consumption = characteristic value of engines 
We follow [13] and model the uncertainties as presented in Table 2. 



variable 


density 


parameter 


V 


Uniform 




F 


Beta 


(7, 2, Ffaim -Fj^iax) 


SFC 




9i ^ 17.23,6(2 = 3.45 



Table 2: Uncertainty modeling 
The probability density function of a beta distribution on [a, b] with shape parameters (a, /3) is 

where B{-, •) is the beta function. Still following [13[, we take the nominal and extremal values of V 
and F as in Table 3. 



variable 


nominal value 


min 


max 


V 


231 


226 


234 


F 


19 


18.7 


19.05 



Table 3: Minimal and maximal values of uncertain variables 
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The uncertainty on the cruise speed V represents a relative difference of arrival time of 8 minutes. 
The airplane manufacturer may wonder whether he has to improve the quality of the engine (SFC) 
or the aerodynamical property of the plane (F). Thus we study the sensitivity of Mf^ei with respect 
to F and SFC and we want to know if Hq : S'^'^^ > or Hi : S'^^^ < . Applying the test 
procedure described previously we can not reject Hq. 

4 Concentration inequalities 

In this section we give concentration inequalities satisfied by the Sobol indices in one dimension (i.e. 
k = 1). 

We define the h function by h{x) = (1 + x) ln(l + x) — x for all x > —1. 

4.1 Concentration inequalities for 5']^q 
We introduce the random variables 

C/± = - ± y)iY,f and J± = {S^^ ± y)Y, ~ F," 

and denote Vjf (respectively Vj7, Vf and Vy) the second moment of the i.i.d. random variable 
(resp. U~, J+ and J~). 

Theorem 1. Let b > and y > 0. Assume that all the random variables Yi and Y^^ belong to [—6, b]. 
Then 

P {S^,ci > Sci + y)< Ah + 2Al2 + 2M3, (10) 
P {S^,ci < Sci ~y)< Mi + 2M2 + 2M5, (11) 



where 



M J ^^UufbuyV\\ ( NV+b\f bu [W^ 

Mr=exp^-^/.(^^-jj ^3=exp|-^.(^^./^ 



1 Wj 



and bu =b'^{l + S^i+y). 

Remark 4.1. One must be cautious since the variables Yi ^Yn are dependent. 

Proof. Since Sq^ and S"" qj are invariant by translation on Y and K", one may assume without loss 
of generality that Y is centered. 

1. Obviously [/^^ and U~ are upper-bounded by b(j, and J~ by bu/b, 

E{U+) = ~yV E(J+) = 
EiU-)^yV E(jr)^0 

and 

= Var(yF") + (5gi + y)'^Vai{Y^) - 2(5'gi ± y)Cov(yF", Y^) + y^V^ 
Vf = ((5Si ± y)' + 1)V - 2(5Si ± y)Cu. 

2. Proof of (10). Using 

{a + b> c} C {a> c/2}U {b> c/2} and {ab > c} C {\a\ > y/c} U {\b\ > y/c} 



11 



one gets 



'(^]^,ci>^a + y) 



7f J2i=i YiY^^ — YnY ^ 



N 



N 



^ ^ (C/+ - E(c/+)) + YnJn > yv] 



N 



< p(^(C/+-E(C/+))]>iV^)+p(5:K.> 



N 



lyV 



Inequality (10) conies directly by applying five times Bennett inequality (see [1] and references 
therein). 

3. Proof of (11). In the same way, one gets 



< p|^f^(_c/-+E(f/-))>iv^^+p|^f]y,>ivy^^ 



Inequality (11) comes directly by applying five times Bennett inequality. 



4.2 Concentration inequalities for T^q 



Now remind Zi = and introduce the random variables Kf = YiY^ — (S'"j ± y) 

Denote (resp. V^) the second moment of the i.i.d. random variable (resp. K~^. 

Theorem 2. Let 6 > and y > 0. Assume that Y G [—6, b]. Then 

P {T^.ci > ^Sl + y)< mi + 2m2l {Sg,+y-l>0}, 

P {n,ci < - y) < mg + 2m4l{sg,+,-i>o}, 

iV(F + C) / 6 ' 



□ 



TOi = exp 



7713 CXp 



62 

"(7 



2 



ni2 = exp ■ 



m4 ~ exp ■ 



262 

N{V + C) 
262 ' 



6 ' ' 



2yV 



V + C\y + 1-S^ 



(12) 
(13) 



ci 
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Proof. Since is invariant by translation on Y and 1'", one may assume without loss of generality 

that Y is centered. 

1. Obvisouly and are upper-bounded by bu, E{K^) ~ ~yV, E{K^) = yV and 



We also have Zi is upper-bounded by b, E(Zi) = and E{Zf) 
2. Proof of (12). One gets if fi-gi + y - 1 > 



2 \ _ y+c„ 



2 



1^1=1 2 l^^ 

= P E (^^'.^ - ii^(^^'^)) + (^^)' (^Si + y - 1) > y^j 

< P ( ^ (A"+ - E{K+)) > iV— ] + P [ ^ > 



2(5Si-Hy-l) 



Inequality (12) comes directly by applying Bennett inequality to the random variables , Zi 
and —Zj. 



3. Proof of (13). One gets since y + 1 - S^^ > 

Pfei<^a-y) = v(^-^J2{-K-+EiK-)) + {ZN)\y + l-S^,)>yvj 



- ^ ( E ("^«" + ^ ) + ^ ( E z^^^ 



2 / y2iy+l-S^,) 

+P ( > N 



2{y+l~S^,)J- 

Inequality (13) comes from Bennett inequality to the random variables , Zi and —Zi. □ 

4.3 Numerical applications 

In this section, we provide numerical illustrations of the concentration inequalities stated in Sections 
4.1 and 4.2. 

The upper bounds appearing in Theorem 1 involve the (a priori) unknown quantities: 

Q^{V,V+.V^,V+,Vi,S^,). 

We denote by pAbove{y, N) and pBelow{y, N) the estimators of the right-hand sides of (10) and (11), 
respectively, obtained by replacing the Q vector by its empirical estimate. 

Similarly, we denote by pAbove' {y, N) and pBelow' {y, N) the estimators of the right-hand sides of 
(12) and (13) when: 

0' = (i-,c,F+,y^,5S,) 

is replaced by its empirical estimate. 

One should note at this point that, on the one hand, the bounds of Theorems 1 and 2 are fully rigorous 
for any N. From a practical point of view, these bounds are not computable, unless the Q (resp. Q') 
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vector is known. On the other hand, pAbove and pBelow (resp. pAhove' and pBelow') are computable 
but are not fuUy justified for finite iV, as they rely on the estimation of Q (resp. Q'). However, as 
pointed out in [4], these bounds are conservative^ hence they are less sensitive to a bad estimation 
than the asymptotic confidence interval given by the CLT 

We again take for / the Ishigami function considered in 3.4.1. 

In this case, it is easy to check that Y G [—6, 5], where: 

6 = 8 + 0.1 X -K^. 

When such a majoration of Y is not possible, h can be put into the Q (or Q') vector and estimator of 
it can be plugged in to obtain p Above and pBelow (or p Above' and pBelow'). 

We also choose u = {1}. 

Figure 3 show, for different values of iV, the plot of pAbove{y, N) and pBelow{y, N) (respectively, 
pAbove' {y, N) and pBelow' (jj, N)) as functions oiy. 




Figure 3: Plots, for N e {1000, 4000, 8000, 12000}, ofpAbove{y, N) (left-top) and pBelow{y, N) (right- 
top). pAbove' {y, N) (left-bottom) and pBelow' {y, N) (right-bottom) for the Ishigami model and for 
u = {l}. 



As expected, the concentration inequalities are more conservative than the asymptotic confidence 
interval. These plots confirm that the concentrates faster than (^j, and the inequality, while 

conservative, is sharp enough for this desirable property of to be reported. We also notice 

that there is a dissimetry in the bounds for above and below deviations, as this is often the case for 
concentration inequalities. Finally, the expected convergence for N — > +oo is observed. 



5 Berry-Esseen Theorems 

In this section we will give a general Berry-Esseen type Theorem for the estimator S" in one 
dimension (i.e. fc = 1). Let $ be the cumulative distribution function of the standard Gaussian 
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distribution. 



5.1 Pinelis' Theorem 



We first recall a general Berry-Esseen type theorem proved in [12]. Let {Vi)i>i a sequence of i.i.d. 
centered random variables in R*^, for some d E K* . Let / some measurable function: M'* R with 
/(O) = and such that: 



3e > 0, 3Afe > s. t. \f{x) - L{x)\ < ^\\x\\^ 
where L := Df{Q) is the Frechet derivative of / at point 0. 



(14) 



Remark 5.1. Remark that condition (14) is satisfied as soon as f is twice continuously differentiable 
in a neighborhood ofO. 

Theorem 5.2 (Corollary 3.7 in [12]). Take any p G (2,3]. Assume (14) holds, 



a := v/E(L(F)2) > 0, 
and (E ([[V^ll^))"'"^^ < oo where \\ ■ \\ denotes the euclidean norm on M''. Then for all z G 



< 



lP/2-l ' 



(15) 



where k above is a generic constant that depends only upon p. 



5.2 Theoretical result for the general case 

For any random variable Z, denote by Z'^ its centered version Z — K{Z). 

Theorem 5.3. Assume that the random variable Y has finite moments up to order 6. Then, for all 
z€R, 

' ' ^ ' '<^. (16) 



a 



N \ 



Here 



is the asymptotic variance of v NS^ 



a' :=Var(l(r'=(yT-5Si(n': 



Proof We define V, = {Yf Cu,Yf (Yf)^ - V) and / : R^ ^ R as f{x,y,z,t) 

T-y'^Vv - ^Ci- Note that /(0, 0, 0, 0) = 0, /(Fjv) = S^ ^i - ^ci and by Remark 5.1, (14) holds. Tl 
result is then a direct application of Theorem 5.2 once cr^ = E (^L{V)'^) will be computed. We have 

df df df df \ , , / 1 

^ ^ (0,0,0,0)= -,0,0, 



^dx^ dy^ dz^ dt J ' ' ' y y ' ' ' y 2 
Using notation in Section 5.1 one gets 

L(x, y, z, t) = 1 [x - 5S,t) and L{V) = 1 {Y^iY^f S^,{Y^f) . 
Straightforward computations lead to the required result. □ 

Then we have a Berry-Essen theorem for Sobol index estimator in a general case (whatever the first 
moment of Y). However, the constant of the bound is hard or even too complex to express explicitely. 
In the next section we present a Berry-Essen theorem with explicit bounds in the centered case but 
with an estimator of S'"j slightly different from q^. 
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5.3 Practical result in the centered case 

In this section we give a Berry-Esseen theorem for the estimator 

— V Y Y^ 



'^N,Cl 



in the centered case and k — 1. Further, let k w 0.42 be the last best constant known in the classical 
Berry-Esseen theorem ([8]). We then have 

Theorem 5.4. Assume that the random variable Y has finite moment up to order 6. Then, for all 
t G R, 



^iS'^,ci'S}^,)<t]-m 



< 



N 



t 



Here 



a' ~ Var ( - (yy" - S^^Y^] 



is the asymptotic variance of v NS^ qj and 



E 



A„ - E(A, 



VVarA„ 



An := CT^V 



t(j 



Y' 



ta 



25^1 ) Var(y2) - 2Cov(yy", y^). 



(17) 



(18) 



Proof. To begin with, we compute the asymptotic variance a of v NS^ : we apply the so- 
called Delta method [17] to = {YiY^,Y^) and ^{u,v) := uv-\ (u G M, v > 0). Then cr^ = 
J*(E(M^))I] J^(E(M^))* (J^ the Jacobian of ^I') and the expression given in (18) follows obviously. 

Now, for i G M, set 



A, 



Obvious algrebraic manipulations lead to 



N 



where, for j = 1, . . . , iV, 



Aat , := (T ^ 



to- 



Y,Y-V-la. + ^V Y 



Now, we have 



C(Ajv) = and Var(Ajv) = 



N 



tVN 



So that denoting by An,- the empirical mean of (Ajvj")j=i,...,7v, we obtain 



|\/iVAAr,. < 0| = |\/7V 



Ajv, -E(AAr) 



N 



< 



Now, to conclude we apply Berry-Esseen theorem (see [8]) and the triangular inequality to obtain 
(16). □ 
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5.4 Numerical applications for the centered case 

We denote by B{t) the right hand side of the Berry- Esseen inequahty (17). It is clear that, for any 
y > 0, we have: 



\-y < s^,ci - s^i < y) > 




and: 



'hv < S^^ci - < y) < 





B I ^y]+B 



N 



(19) 



(20) 



Hence, the actual confidence level of the asymptotic confidence interval for S*" using S]^ q is greater 
than the theoretical level (first term of the sum above), minus a correction term given by the Berry- 
Esseen theorem (second term). The upper bound given by (20) may also be of practical interest: an 
overly conservative (overconfident) interval is not always desirable, as a more precise interval with 
accurate level may exist. 

As in the previous applicational section 4.3, the lower bound of the asymptotic confidence interval 
level involve unkown quantities (moments of An, Y, YY^) that have to be estimated. We designate 
by L{y,N) (resp. U{y,N)) the estimator of the right hand side of (19) (resp. (20)) when all unkown 
quantities are empirically estimated. 

We take as output model the Ishigami function defined at Section 3.4.1, recentered by its true mean 
7/2 : 



Y = f{Xi,X2, X3) = sinXi + 7 sin^ X2 + O-IX^ sinXi 
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Note that the true mean could also be replaced by an estimate of the mean. For y, we choose 
y = 1.96-^^, where is an empirical estimate of cr^, so as to compute (estimators of ) upper and 
lower bounds of the actual level of the 95%-level confidence interval. 

We present the numerical results, as functions of N, and for u = {1} in Figure 4; for u 
u = {3}, the results were very similar. 



{2} or 



1.3 



1.2 



7T 1.1 



0.6 




100 



Minoration of CI level 
Majoration of CI level 
Asymptotic level 



^ 



1000 



10000 100000 
N 



1e+06 



1e+07 



Figure 4: Plots of L{N) (minoration of CI level)) and U{N) (majoration of CI level) for u = {1} and 
different values of N. 



As expected, the actual confidence level is estimated under the "target" level of the confidence interval 
(0.95). As +00, our bound converges (quite slowly) to 0.95. Nevertheless, the Berry-Esseen 
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bound we have presented quickly attains confidence levels which are very close to the asymptotic level, 
and it can be used so as to provide a certification, at finite sample size, of the level of the asymptotic 
confidence interval. 
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